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ABSTRACT 

Recent observations have shown that at least some close-in exoplanets main- 
tain eccentric orbits despite tidal circularization timescales that are typically 
much shorter than stellar ages. We explore gravitational interactions with a 
more distant planetary companion as a possible cause of these unexpected non- 
zero eccentricities. For simplicity, we focus on the evolution of a planar two-planet 
system subject to slow eccentricity damping and provide an intuitive interpreta- 
tion of the resulting long-term orbital evolution. We show that dissipation shifts 
the two normal eigenmode frequencies and eccentricity ratios of the standard sec- 
ular theory slightly, and that each mode decays at its own rate. Tidal damping 
of the eccentricities drives orbits to transition relatively quickly between periods 
of pericenter circulation and libration, and the planetary system settles into a 
locked state in which the pericenters are nearly aligned or nearly antialigned. 

Once in the locked state, the eccentricities of the two orbits decrease very 
slowly due to tides rather than at the much more rapid single-planet rate, and 
thus eccentric orbits, even for close-in planets, can often survive much longer 
than the age of the system. Assuming that an observed close-in planet on an 
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elliptical orbit is apsidally-locked to a more distant, and perhaps unseen compan- 
ion, we provide a constraint on the mass, semi-major axis, and eccentricity of the 
companion. We find that the observed two-planet system HAT-P-13 is in just 
such an apsidally-locked state, with parameters that obey our constraint well. 
We also survey close-in single planets, some with and some without an indication 
of an outer companion. None of the dozen systems that we investigate provide 
compelling evidence for unseen companions. Instead we suspect that i) orbits are 
in fact circular, ii) tidal damping rates are much slower than we have assumed, 
or iii) a recent event has excited these eccentricities. Our method should prove 
useful for interpreting the results of both current and future planet searches. 

Subject headings: Planetary Systems, Planets and satellites: dynamical evolution 
and stability, Planets and satellites: general 



Introduction 



In the past decade, many mechanisms have be en proposed to ex plain a wide range of ec 



Namounil ( 120071 ) and references there-in) 



centricities seen among exoplanets (e.g., a review by[ 
including planet-planet scattering, planet-disk interaction, mean-motion resonance passage, 
and Kozai resonances. All these mechanisms can excite orbital eccentricities effectively. 
While the orbits of long-period planets could stay eccentric for billions of years, those of 
close-in planets are likely to be circularized within the stellar ages due to tidal interactions 
between the stars and the planets. 

In Fig. m we plot the eccentricities of all 222 known (as of February, 2012) close-in plan- 
ets with orbital periods less than 20 days. The average eccentricity for these planets (~ 0.06) 
is significantly smaller than that of all the confirmed planets (~ 0.18) which indicates that 
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Fig. 1. — Eccentricities of close-in planets with orbital periods less than 20 days. The blue 
circles and orange triangles correspond to single- and multiple-planet systems, respectively. 
Error bars are plotted for both period and eccentricity. However, errors in period are not ap- 
parent because they are very small. Orbital data, courtesy of the Exoplanet Orbit Database 



( http:/ / exoplanets.org/ ) . 



the tidal interaction is an effective eccentricity damping mechanism, and that circularization 
timescales are typically shorter than stellar ages (about 1-10 Gyr). Nevertheless, nearly half 
of these close-in planets have non-zero orbital eccentricities. Recent studies have shown 



however, that t 


le or 


Zakamska et al. 


2011 



Shen fc Turner 



2008 



Pont et al. 



201ll ). so some measured eccentric orbits for close-in plan- 
ets may, in fact, be circular. But it is unlikely that all close-in planets have perfectly circular 
orbits and thus at least some of the non-zero eccentricities require a physical explanation. 
The main possibilities include i) systems have only recently attained their orbital configu- 
rations and eccentricities are damping quickly, ii) planetary tidal quality factors are much 
larger than those of the giant planets in our Solar System and so eccentricities damp slowly. 
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and iii) eccentricity excitation due to an exterior planet slows the orbital circularization. In 
this paper, we explore the third option and investigate the gravitational interactions between 
a close-in planet and a more distant companion. 

In systems with more than one planet, the most significant orbit-orbit interactions are 
often mean-motion resonances, whi ch have been studied in deta il for the satellite systems 
of the giant planets (see reviews by 



Greenberg 



1977 



Peak 



planetary 



2003 



Lee 



systems have also received increased attention (e.g., 



2004; 



Ketchum et al. 



2013 



Batygin fc Morbidelh 



1986). Resonance s in extrasolar 



Chiang 



2003 



Beauge et al. 



20131 ). Mean- motion resonance 



passages during planetary migration can be effective in exciting orbital eccentricities. These 
resonance passages can be divergent, in which case impulsive changes to the o rbital ele 



ment s result, or convergent, in which case trapping into resonance usually occurs (IHamilton 



1994 



Zhang &: Hamilton 



2007 



20081 ). Although resonance capture typically leads to ex- 
cited eccentricities, we do not consider this process further here. Instead, we focus on the 
more prosaic secular interactions which are capable of maintaining orbital eccentricities for 
a greater variety of orbital configurations. 

Secular p erturbations have be en studied for centuries in the context of the Solar Sys- 
tem (see, e.g. . 



Butler et al. 



2002 



2009 



Brouwer et al. 



19501 ). With the discovery of extrasolar multi-p 



19991) ■ apphcations 



Barnes fc Greenberg 



2006b 



br secular theory haye expanded rapid 



Greenberg Sz Van Laerhoven 



Adarns fc Laughlin 



2011 



Laskar et al. 



2006a 



Mardling 



y (e.g . 



2007 



anet systems 



Wu fc Goldreich 



Batvgin et al. 



Wu fc Goldreich! fl2002h were 



the first to use secular interactions to explain the non-zero eccentricity of a hot Jupiter. They 
showed that a close-in planet with a compani on could sustain a substantial orbital eccentric- 



ity even though tidal damping was efficient. 



Mardling] ( 120071 ) developed an octopole-order 



secular theor y to study orbita l evolu tion of close-in planets with companions and confirmed 



the results of 



Wu fc Goldreich! (120021 ). If secular interactions play an important role in delay- 
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ing eccentricity damping, we may expect differences in the eccentricity distributions between 
single- and multiple-planet systems. However, as seen in Fig. [H there is no obvious difference 
between these two groups of planets. Although future observations may reveal such a differ- 
ence, it is also plausible that many single close-in planets are accompanied by unobserved 
companions that help maintain their eccentricities against tidal dissipation. 

In this paper, we revisit the problem of secular interactions with a distant companion 
in maintaining the eccentricities of close-in planets. Our goals are to develop an intuitive 
interpretation of the secular theory of a two-planet system and to test our model against 
observed planetary systems. In the next section, we present our model for secular orbital 
evolution during tidal dissipation, starting with a review of secular orbital interactions in 
a stable non-dissipative system consisting of a star and two planets. We then add tidal 
dissipation of eccentricities, and solve the coupled system to investigate how eccentricity 
damping affects the apsidal state of the two orbits. Finally, we add the effect of General 
Relativity (GR) which can be significant for close-in planets. In Section [31 we apply our 
model to extrasolar planetary systems, illustrating how our model may help to guide planet 
searches. Finally, we discuss and summarize our work in Section HJ 



2. Model 



2.1. Stable Non-dissipative Systems with Two Planets 



The secular solution of a stable non-d issipative 2-planet system has been known for 



centuries and is discussed in great detail in 



Murray fc Dermott 



( I1999I ). In this section, we 



will develop a graphical interpretation of the solution that will help us understand the more 
complicated systems studied later in this paper. After averaging the disturbing functions 
of each planet on the other, TZj, over both planets' orbital periods, Lagrange's planetary 
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equat ions can be linearized for small eccentricities and inclinations to (IMurray fc Dermott 



1999 



, §7.1) 



dj = 0; 



1 dUj . _ 1 dTZj 

1 dn. 

I 



Here rij, Oj, e^, Ij, flj, and ccj are the mean motion, semi-major axis, eccentricity, inclination, 
longitude of ascending node, and argument of the pericenter of the j^^ planetary orbit, 
respectively. 

Planetary semi-major axes remain constant and hence no long-term energy transfer 
between the orbits occurs because we have averaged a conservative perturbation force over 
each planet's orbital period. Note that, due to our assumption of small eccentricities and 
inclinations, the evolution equations of (ej, Wj) and {Ij, Qj) are completely decoupled and 
can be analyzed separately. In this paper we neglect the inclination and node pair since they 
are less easily observable for extrasolar planets. Because the two sets of equations take the 
same form, however, our eccentricity results below can be easily applied to secular coupling 
of vertical motions. 

The dis turbing function 7^ is defi ned for a planet as the non-Keplerian potential at 



its location ( Murray fc Dermott 



19991 . §6). For two-planet systems without any external 



perturbations, it is simply the gravitational potential due to the other planet. Here we 
consider a simple planetary system consisting of a c entral star and two planet s in co-planar 
orbits. In terms of osculating ele ments (see, e.g., 



second order in small eccentricities, 
disturbing functions are: 



Murrav fc Dermott 



Murray fc Dermott 



1999 



^2.9) and to 



19991 ) show that the orbit-averaged 



TZi = nia\ aq 



-el — (3 6162 cos(ct7i — W2) 



(2) 
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7^2 = ^^202 



1 



(3 6162 cosfroi — -072) 



(3) 



where the subscript "1" refers to the inner planet and "2" to the outer one. We define the 
mass ratio between the two planets q = 7712/7711, and the semi-major ratio of the two orbits 
a = ai/a2- The remaining parameters are defined as: 

M 1 mi 2(1) 



3/2^ 



and a = -7ii — a 6.;,!, (a) 
4 ■^1'^ 



w here m^. is the s te^ 



( iMurray fc Dermott 



lar m ass and 63^2(0) and 63^2 («) are two of the Laplace coefficients 



19991 ). The parameter a has units of frequency, which we will soon see 



characterizes the secular precession rates. Both /3 and a decrease with increasing planetary 

'3/2 (") ~ 3a and h^^j^{ 

Following 



separation and, for small a, reduce to h^}}ja) ~ 3a and h^?}ja) ^ 15a^/4. 



Brouwer et al. 



( 1l950l ). we transform a (ej, Wj) pair into a complex Poincare 



canonical variable hj with the mapping: 



J/' 



(4) 



where i = Substituting Eqs. (I2]|3]) into Eqs. ([T]) and rewriting in terms of hj yields a set 

of linear homogeneous ordinary differential equations similar to those for a double pendulum 
system: 



hj = i^Ajkhk, 

k=l 



(5) 



where the coefficient matrix 




A 



Two orthogonal special solutions of Eq. (|5]), or the secular eigenmodes of the system, 
are given by 













^ v± ) 



exp(i^^t). 
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Aligned Mode 



Anti-Aligned Mode 




Fig. 2. — Aligned and anti-aligned secular modes. The planets (solid dots) follow nearly 
elliptical orbits about the central star. Arrows point from the star to orbital pericenters. 

where the eigen-frequencies and eigenvector parameters rj^ can be obtained from the 
matrix A: 



9± 



a 
2 



q + y/aT \ {q- Va)^ + 4:qy/a(3'^ 



q-^±^{q- ^/^Y + 4gv/^/j^ 
2q(3 



(6) 
(7) 



We use the superscript "s" to indicate that those parameters are for a "static", or non- 
dissipative, system. Notice that we have chosen our eigenmode components hi± and /i2± to 
depend only on the fixed constants a and q through Eqs. [6] and [71 and not on the initial 
eccentricities and pericenter angles. 

The physical meaning of the two modes can be elucidated by transforming the solution 
of hj back to the orbital elements (e, w) with Eq. (jl]). If the system is fully in either the 
"+" or the "-" mode, we find 



■wi± = zu2± = 
(e2/ei)± = \v±\, 
cos(Acc7±) = ?7^/|?7^| 



(9) 
(10) 



where Aw = W2 — ^i is the difference between the two pericenter angles. In either mode, the 
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two orbits precess together at the rate g± (Eq. [8]), and their eccentricities keep a fixed ratio 
(Eq. [9]). Furthermore, Eq. ([7]) shows that 77^ > while rjt < 0. Thus, Eq. ffTOj) states that 
the pericenters of the two orbits are always aligned in the "+" mode (cos(A?z7) = 1), and 
anti-aligned in the "-" mode (cos(Ati7) = —1). In fact, in an eigenmode, the system behaves 
as a rigid body with the shapes and relative orientation of the elliptical orbits remaining 
fixed (Fig. [2]). The frequency of the anti-aligned mode is faster {g_ > g^) due to closer 
approaches between the planetary orbits and thus stronger perturbations (Fig. [2]). 

It is instructive to consider the small q limit which corresponds to a tiny outer mass. In 
this case, Eqs.[6]and[7]simplify to {g^ = crg(l— 77^ = (3) and {gt_ = ay/a, 77! = —y/a/qP). 
For a true outer test particle g — )■ and only the first mode is possible. This aligned mode 
is stationary ((7+ = 0) and, since 77^ < 1, the planet's eccentricity exceeds that of the test 
particle. Similarly, with an inner test particle only the aligned mode survives, it is stationary, 
and the massive planet again has the higher eccentricity. 

Eq. ([9]) gives the ratio between the eccentricities in a secular mode. In Fig. |3l we plot 
77^ as a function of y/a/q for different (3 values. Although rare in real systems, q = y/a (or 
mlai = m\a2) makes an interesting special case. When this condition is met, the planets 
have similar angular momenta and the inner and outer planet orbital precession rates (the 
diagonal terms of the matrix A) are equal. In addition 7/j_ = ±1 (Eq. [7]), and therefore, the 
inner and outer orbits have the same eccentricity (Eq. [H]). When q > i/a, the inner planet 
precesses fastest, has a lower eccentricity in the aligned mode and a higher eccentricity in 
the anti-aligned mode. The opposite is true for q < a/o {mfai > 7722(22) • 

In general, a system is in a mixed state composed of a linear combination of the two 
modes: 

hj = e+ exp(i<y9+) hj-^ + e_ exp(i<y9_) hj-. (11) 
Here the mode amplitudes e± and phases ip±. are determined by the initial orbits. Plotting 
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2 

»+ 1 


s=- -1 

-2 

0.0 0.5 1.0 1.5 2.0 

^/a|q 

Fig. 3. — Eccentricity ratios versus \/a/q in the secular modes as given by Eq. (jTj). Different 
curves represent different /3 values. The dashed line (g = ^/a or m^ai = 772202) is a turning 
point; these orbits have ei = 62 in both modes. 

the solution Eq. (fTTj) on the complex plane yields a phase plot of ecoscc versus e since as 
shown in Fig. UK- The ei and 62 vectors in the plot represent (e, vj) pairs for the two orbits 
at a given time. The length of a vector is the instantaneous eccentricity, and its polar angle 
is the instantaneous longitude of pericenter. Each eccentricity vector is a vector sum of an 
aligned component (e^ and ?7^e+ are parallel) and an anti-aligned component (e„ and ?7ie_ 
are antiparallel). The lengths of all of these vectors are determined by initial conditions. 
In a pure aligned eigenmode, e_ = |e_| = and the ei and 62 vectors are parallel. As 
time progresses, these eccentricity vectors rotate together at rate while maintaining their 
lengths. This solution corresponds to the aligned orbits in Fig. [2j For the anti-aligned 
eigenmode, e+ = |e_|_| = so that ei and e2 are anti-parallel and rotate together at rate g'L. 
This state is depicted by the pair of anti-aligned orbits in Fig. |2l In the most general system, 
both motions occur simultaneously: the parallel eccentricity vectors (e+ and ?7+e+) rotate 
at rate f?^, while the anti-parallel vectors (e_ and ?7ie_) rotate at rate g'L- The resulting 
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Fig. 4. — General solution for a two-planet system, a) The solution (Eq. [TT]) on a phase 
plot. The arrows represent the eccentricity vector eexp{iuj). The total eccentricity of each 
orbit (ei or 62) is the magnitude of the vector sum of aligned (+) and anti-aligned (— ) 
components. The two aligned components, e_|_ and ?7+e+, rotate (precess) at rate g^, while 
the two anti-aligned components (e_ and rjte^) rotate at rate gt- b) Secular evolution 
of orbital eccentricities from an N-Body simulation. Plot shows the eccentricities of two 
planetary orbits in a computer-simulated system consisting of a 1 Solar-mass star, a 1 Jupiter- 
mass "hot- Jupiter" at 0.05 AU, and a 0.8 Jupiter-mass companion at 0.2 AU. The simulation 
shows an oscillation period of ~ 1673 years, in good agreement with the prediction of the 
secular model: 2n/{g_ — g^) = 1670 years. 
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lengths of the eccentricities ei and 62, thus vary periodically as illustrated in Fig. |3)d. The 
maximum value of e\ occurs when e_|_ is parallel to e_. At the same time, however, ?7^e_|_ and 
77! e_ are anti-parallel to each other, leading to a minimum value for 62- The simultaneous 
maximum for ci and minimum for 62 seen in Fig. |Dd is a general result guaranteed by angular 
momentum conservation. The two eccentricities, in mathematical form are: 



ei = -v/ + e?. + 2e+e_ cos(5fi - g%)t, 
62 = \j (e+?7+)^ + {e-rjLy + 2e+e_?7+?7i cos{gt — g'^)t 



( Murray fc Dermott 



I999I ). Both eccentricities oscillate at the same frequency (^i — g^) as 



shown in Fig. |3)d. 



2.2. Secular Modes with Eccentricity Damping 

Since a planet close to its host star experiences a drag force due to planetary tides raised 
by the star, we seek a way to include tides into the mathematical formalism of the previous 
section. For small eccentricities for which the secular solution Eq. (fTT]) is valid, tidal changes 
in a are usually negligible compared to the damping in e. Since the damping timescales 
for close-in planets (in the order of 10^ yrs) are much longer than the secular timescales 
(typically ~ 10"^ yrs), we treat the damping effect as a small perturbation to the secular 
solution. 



Stars and planets raise tides on each other, which in turn perturb the planet's orbit. For 
very close-in planets, both stellar tides raised by the planet, and planetary tides raised by the 
star act to decrease the planet's orbital period and eccentricity. Typically, tidal dissipation 
in the star leads to orbital decay and circularization on a timescale much longer than the 

dam p an orbital eccentricity 



rather quickly (e.g. 



Burns 



1977 



while that 


in t 


Goldreich 


1963 



Rasio et al 



19961 ). Thus, we assume that 
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the orbital circularization is largely dom inated by planetary tide s, and approximate the 



eccentricity damping rate as follows (e.g. 



Murray fc Dermott 



19991 ) : 



e 63 1 m^. / R. 



n . 



(12) 



where Q'^ = L'^Q^jk^ is the modified tidal quality factor, is the tidal quality factor, and 
k2 is the Love number of degree 2. 

The addition of the constant tidal damping (Eq. [T^ adds an extra term to the eccen- 
tricity equation in Eq. ([T]), which now reads 

. _ _ 1 dTZj _ 

njttjej owj 

with A J being the eccentricity damping rate. The coefficient matrix of Eq. ([5]) is now 



A = a 



where the dimensionless = \j/ a <^ 1 parameterizes the damping strength. 



(13) 



Eccentricity damping causes the eigen-frequencies of matrix A to have both real and 
imaginary parts. As in other dynamical systems, the real parts {g±) of the eigen-frequencies 
still represent the precession rates of the secular modes, while the imaginary parts (7±) 
indicate that the amplitudes of the modes change over time. This becomes clearer if we 
rewrite the two orthogonal special solutions of the system 

(14) 









V h2± j 







exp(-7±t) exp(i5(±t), 

where rj± is the new eccentricity ratio for each mode. 

For ^ 1, we solve the new matrix for the complex frequencies and find 

9± = 9i± 



7± 



[(g- v^)2 + 4gy^/32]3/2 
y/a-q 



Ai + Aa ± 



v/(g- v^)2 + 4gv^/32 



^(^1-6)', 

(Ai - A2) 



(15) 
(16) 
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Eccentricity damping increases the precession rate of the ahgned mode and decreases that 
of the anti-ahgned mode, but only by a n extremely s r nall a mount (of order ^j). These 



tiny frequency changes were neglected by 



Laskar et al. 



(120 121 ) . but our damping rates are 



in perfect agreement. The tiny frequency changes are due to the slightly different orbital 
configurations in the secular modes as we shall describe below. The new eigen- vectors of the 
matrix A are: 



V± = V± 



1 ±t- 



where we have neglected terms of second- and higher-order power of C,j- 



(17) 



In general, the ?7±'s are complex with small imaginary components. If we ignore the 
imaginary parts for the moment, then ?7± are real and Eq. f|T^ shows that the pericenter 
angles of the two orbits are the same (positive ?7+, aligned mode) or 180° apart (negative 
?7_, anti-aligned mode). For complex ri±, however, the two orbits are not exactly aligned or 
anti-aligned any longer. Instead, Aw± shifts from 0° and 180° by a small angle 



e = tan ^ 



V(g- V^)' + 4gv^/3V V(g- v^)2 + 4gv^/32' 



In the "aligned" mode, the new pericenter difference is Aw+ = e so that the inner 
exoplanet's pericenter slightly lags that of the outer exoplanet. The lag is maximized for 
q = i/a? the case with equal eccentricities and equal precession rates for the two planets. 
Because of this mis- alignment, the minimum distance between the two orbits is slightly less 
than that in the undamped case (see Fig. [2]). This causes the average interaction between the 
two orbits to be stronger, leading to an increase of the precession frequency as indicated by 
Eq. ( ITSl) . Similarly, Aro_ = 180° — e in the "anti-aligned" mode; the slight rotation results in 
a weaker average interaction and a slightly slower mode-precession rate. The deviation angle 
e is tiny, and the eccentricity ratios in the two modes, |?7±|, are nearly the same as \ri^\- Thus, 
we continue to use "aligned" and "anti-aligned" to refer to the two modes. Nevertheless, 
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the small deviation angle is physically important because it enables damping of the outer 
planet's eccentricity. 

In addition to the slight mis- alignment, each mode amplitude also damps at the rate 
given by Eq. f ITB]) . If only planetary tides contribute to eccentricity damping, Eq. f[T^ shows 
that the damping rate decreases rapidly with the planet's semi-major axis (A oc a~^'^). In 
the absence of secular interactions between the planets, the outer orbit is hardly affected. 
With this interaction, however, the damping applied to the eccentricity of the inner orbit 
is partially transmitted to the outer planet, causing a decrease of its eccentricity as well. 
The damping rates of the two modes are different, unless q = ^Ja. An interesting result 
from Eq. (fT6|) is that the sum of the two mode-damping rates is equal to the sum of the two 
individual eccentricity damping rates: 

7+ + 7- = Ai + A2. 

The physical interpretation of this expression is that secular interactions between the 
planets simply act to redistribute where the damping occurs. 

In Fig. [5l we compare our analytical results with numerical integration of both the 
secular equations and the direct N-body equations, with an artificial eccentricity damping 
added only to the inner planet in all cases. The eccentricity evolution curves of the inner 
planet are plotted for two different cases, and e-folding rates are measured and labeled for 
all curves. Fig. [5^ illustrates a system with q > ^/a. The top panel shows results from 
secular equations, and the bottom panel shows the corresponding N-body simulations. Each 
panel plots three curves: i) the single planet case, in which the eccentricity of the inner 
planet damps at rate Ai, ii) a two-planet aligned mode with damping rate 7+, and iii) a 
two-planet anti-aligned mode (damping rate 7_). For q > ^/a, Fig. |5ti, eccentricities damp 
much faster in the anti-aligned mode than in the aligned mode, as predicted by Eq. (fT6|) . A 
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q = 0.8, a = 0.25 q = 0.3, a = 0.25 




2xl0« 4xl0e 6xlO« 8xl0» 10' 2x106 4x106 6x106 BxlO» 10' 



t (yr) t (yr) 

(a) (b) 

Fig. 5. — Eccentricity damping of systems in different secular modes found by integration 
of the secular equations (top panels) and direct N-body simulations (bottom panels). The 
plots show the eccentricity evolution of a "hot- Jupiter" (1 Jupiter-mass planet at 0.05 AU 
from a 1 Solar-mass star) with a companion; 7+ and 7_ represent the mode damping rates, 
which are e-folding times measured directly from the curves. Also plotted is a single "hot- 
Jupiter" subject to an artificial eccentricity damping with a rate Ai = 7.90 x 10~^ yr~^- (a) 
A 0.8 Jupiter-mass companion is located at 0.2 AU (g > y/a), with predicted mode damping 
rates (Eq. ([16]): 7+ = 1.5499 x 10"^ yt~\ 7_ = 6.3501 x 10^^ yr^^ (b) A 0.3 Jupiter-mass 
companion is located at 0.2 AU (g < y/a), with predicted 7+ = 6.4779 x 10""^ yr~^ and 
7_ = 1.4221 X 10-^ yr-i. 
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comparison between the top and bottom panel shows close agreement (within 2%) between 
full-scale N-body simulation and integration of the approximate secular equations. Damping 
rates predicted by Eq. f|T6|) match the observed secular decay rates almost perfectly. Fig. [5h> 
shows a system with q < ^/a, for which the aligned mode damps faster than the anti-aligned 
mode. The faster damping rate of the aligned mode in the N-body simulation is still within 
2% of the prediction, but that of the slow anti-aligned mode, however, is ~ 10% off. This 
discrepancy may be due to unmodeled perturbations to the inner body's semi-major axis. 

The different damping rates for the two modes are particularly interesting, especially 
for well-separated nearly-decoupled orbits for which a and /3 are small. In this case, the 
Aq^/aP'^ term under the square root of Eq. (fT6|) is much smaller than the other term. As a 
result, if the eccentricity damping on one orbit is much faster than on the other (Ai ^ A2), 
as in the case of tides, one mode damps rapidly at nearly the single-planet tidal rate Ai. The 
system evolves quickly into a single mode which decays substantially more slowly. 

Due to the different damping rates for the two modes, a system will evolve into a single 
mode even if it starts in a mixed state. Fig. |6] shows the eccentricity and apsidal angle 
evolution of the systems depicted in Fig. [5l but with initial conditions that leads to mixed 
states. Fig. [6^ shows the case of g > y/a. Before 4 Myrs, the system is in a mixed state, 
so both eccentricities, as well as their ratio, oscillate (cf. Fig. H]). As the short-lived anti- 
aligned mode damps away, the orbits begin to librate around Aw ~ 0°, the two eccentricities 
oscillate less and less, and in the end, the eccentricity ratio settles to the aligned mode ratio 
|?7_i_| predicted by Eq. flTTI) . Fig. shows the corresponding plots for q < ^/a. The orbital 
elements undergo similar evolution, except that the aligned mode damps quickly and the 
system ends up in the anti-aligned mode. 
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(j = 0.8,a=0.25 = 0.3, a=0.25 




time (Myr) time (Myr) 

(a) (b) 

Fig. 6. — Secular evolution of the same systems shown in Fig. |5l but with different initial 
conditions so that the systems begin in mixed states, (a) For q > ^/a, the anti-aligned 
mode damps quickly at rate 7_ from Fig. [5^, and the system evolves to the aligned mode 
(Atu ^ 0°). b) For q < ^/a, the aligned mode damps more rapidly and the system evolves 
to the anti-aligned mode {Azu ~ 180°). 



- 19 - 



2.3. Apsidal Circulation and Libration 

In Fig. [6l the apsidal motion of the two orbits changes from libration of Aro about 0° or 
180° to circulation of Au; through a full 360°, and to libration again during the eccentricity 
damping. In order to understand what determines the apsidal state of the orbits, we plot the 
aligned and anti-aligned components of the eccentricities for Fig. [6^ in the top two panels of 
Fig. [71 Recall that the total eccentricities of the orbits at any time can be obtained from the 
components as illustrated in Fig. SJ The lower panels in the figure show the phase diagrams 
of both orbits (ecos(AG7) versus esin(Aro)) at five different time indicated in the top two 
panels. The two orbits move along the phase curves, which themselves change slowly over 
time. Two critical instants, labeled s_ and s^, are circulation-libration separatrices, which 
represent the transitions of the apsidal state from anti-aligned libration to circulation (s_), 
and from circulation to aligned libration These two points divide the evolution curves 

into three regions. 

In region 1 (t < 10^ yrs), the anti-aligned components are stronger than the aligned 
ones for both orbits (e_ > e+ and |?7!_|e„ > ?7^e+). In the phase plots, both the ei and 
the 62 curves are closed and stay on the left side of the vertical esin(Acc7) axis, indicating 
the libration of Aw about 180°. With the decrease in the amplitudes of all components, 
especially the faster damping of the anti-aligned ones, the curves move closer toward the 
origin, resulting in an increased libration width of Azu. 

The anti-aligned separatrix s_ crossing occurs at t = 10^ yrs, when the two components 
for the outer orbit are equal (|?7!.|e_ = ?7+e+) and 62 may drop to zero, resulting in a phase 
curve for the orbit that is tangent to the vertical axis at the origin (s_ curve on the inner 
phase plot of Fig. [7]). The phase curve for ei at s_ is a half-oval whose straight edge includes 
the origin. Accordingly, when 62 drops to zero. Aw jumps from 90° to —90° for the largest 
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Fig. 7. — Evolution of the apsidal states during eccentricity damping. The top two panels 
show the time evolution of the eccentricity components for the system in Fig. [6^ for which 
mfai < m2a2 and the anti-aligned mode damps fastest. The inner orbit has an aligned 
component e+ and an anti-aligned one e_, while ?7+e+ and |?7!.|e_ are the components for the 
outer orbit. Two circulation-libration separatrices (s_ and divide the evolution curves 
into three parts: anti-aligned libration (region I), circulation (region II), and aligned libration 
(region III). The bottom panels show phase diagrams of the inner and outer orbits on the 
complex eexp{iAw) plane at the corresponding points indicated in the top panels. The 
shape of the phase curves depends on the relative strength of the two components for each 
orbit. Banana shapes result from the large difference between the two components: e_ <^ e+ 
at (a) and |?7f_|e_ <^ f]+^+ (c)- 
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possible full libration amplitude of 180°. 

As the system moves past s_, the phase curves for both orbits enclose the origin and 
circulation results. The circulation region (II) is located between the two separatrices (i.e., 
when 10^ yrs < t < 3.95 x 10^ yrs), where the anti-aligned component of the inner orbit is 
stronger than the aligned one (e_ > e+), while it is weaker for the outer orbit (|?7l|e__ < 
?7^e+). With the continuous fast damping of the anti-aligned mode, the system crosses the 
aligned separatrix 8+ at t = 3.95 x 10^ yrs, when the two components for the inner orbit are 
equal (e_ = e+). The two separatrices occur at those times when each phase curve in Fig. [7] 
touches the origin. 

After s+, both phase curves are to the right of the vertical axis, indicating libration 
about the aligned mode (region III). The two anti-aligned components are both significantly 
damped and the system now has both e_ < e_|_ and |?7!_|e_ < ?7^e+. 

The geometry of the orbits can also be illustrated with a component diagram similar to 
Fig. m but in a frame rotating at the same rate as the anti-aligned mode (Fig. [8]). Note that 
in this rotating frame, the aligned component vectors always rotate clockwise because their 
precessions are slower than those of the anti-aligned ones. Evolution in these coordinates 
can be visualized as circles whose radii and horizontal distances from the origin shrink at the 
different rates 7+ and 7_. Since the anti-aligned components initially dominate the aligned 
ones (region I), the ei vector stays on the right side of the vertical axis, and the e2 vector 
on the left side in Fig. |H1 When the aligned components are parallel to the vertical axis, 
the angle between ei and 62, |Ati7| > 90°, is at minimum; thus the orbits librate about 
Agj = 180° (anti-aligned libration). When the 62 circle moves to enclose the origin and the 
ei circle is still confined in the first and fourth quadrants, the system reaches circulation 
region II. This geometry enables Aw to cycle through a full 360° (Fig. |8]). Finally, when 
the anti-aligned components are sufficiently damped and both circles contain the origin, the 
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Fig. 8. — Eccentricity component diagrams for different regions in Fig. [71 These diagrams 
are similar to Fig. HJ but now shown in a frame rotating at the same rate as the anti-ahgned 
mode so that the horizontal vectors are e_ and ?7_e_. Here the anti-aligned mode damps 
faster than the aligned mode (7_ > 7+) so that the circles move horizontally toward the 
origin faster than their radii shrink. The system starts with Aro ~ 0° in region I and evolves 
to Ati7 ^ 180° in region III. 
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system goes to aligned mode libration (region III). Now the maximum value of Atu < 90° 
occurs when the two aligned components are parallel to the vertical axis, and the orbits 
librate about Atu = 0°. 

Fig. [H] shows the case of Fig. Eb, where the two aligned components are initially stronger 
and the system starts in the aligned libration region III. The two orbits evolve to cross the 
aligned separatrix s+ into the circulation region II, and then pass the anti-aligned separatrix 
S- to reach the final anti-aligned libration region I. The equivalent of Fig. [8] for this system 
would show two circles that initially encompass the origin; here the radii of the circles would 
shrink faster than the distances of their centers from the origin. 

In conclusion, the apsidal state of a two-planet secular system depends on the sign of 
the simple product 

g = (e+-e_)(r/;e+ + <e_) > 0. 

Libration occurs when the same mode components are stronger for both orbits (Q > 0), and 
circulation occurs when one mode is stronger for the inner orbit, but weaker for the outer 
on e (Q < 0). This result is in fu ll agreement with a slightly more complicated formula given 



by 



Barnes fc Greenberd (l2006bl ). 



Eccentricity damping is effective in changing the apsidal state of the orbits because 
the two modes damp at different rates. In fact, eccentricity excitation is equally capable of 
moving the two orbits across libration-circulation separatrices. This can be easily visualized 
by running the plots in Fig. [7] and M backwards in time. For systems with m^ai < 772202 
(see Fig. [7]), eccentricity excitation would eventually bring the orbits into anti-aligned libra- 
tion (region I), while eccentricity damping brings them into aligned libration (region III). 
The opposite is true for systems with mfai > m\a2 (Fig. |9]). All mechanisms that change 
eccentricities slowly cause planetary systems to move toward apsidal libration. 




Fig. 9. — Evolution of the apsidal state during eccentricity damping. Similar to Fig. [71 but 
using data from Fig. |6b for which mfai > m\a2- The aligned mode damps fastest and the 
system moves from aligned libration to circulation, and finally to anti-aligned libration. 
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2.4. General Relativistic Correction 

For close-in exoplanets, the general relativistic (GR) effects are important. In particular, 
it is well-known that the post-Newtonian potential pro duces an apsidal precession of the orbit 



which, to lowest order in eccentricity, is (jPanby 



19881) 



2^3 



3a^n 

^GR — 



where c is the speed of light. We define a dimensionless quantity n to measure the relativistic 
effect: 

K = , 

a 

which is the ratio of the relativistic precession to the characteristic secular precession. 
General relativity adds an extra term to Eq. ([T]), which now reads 

Since rij*^^ is independent of e for small eccentricities, the extra terms do not change the 
form of Eq. ([5]), and so all discussion of the general secular modes still holds. In particular, 
the system still has aligned and anti-aligned modes and the two modes damp separately. 
The mode frequencies, damping rates, and eccentricity ratios, however, need to be revised. 
Now the diagonal terms of the coefficient matrix Aj^ should be adjusted to 



A = a 



q + K- + i^i -qf3 

which gives the new mode frequencies and eccentricity ratios: 

= ^ o- 1 (g + k) + ^/a{l + a\) :^^[q + k- ^/a{l + a'^n)]'^ + 4 qy/a(3'^ 

^ _ q + K- ^/a(l + o^k) ± yj\q + k - ^/a{l + a^K)Y + Aq^l3'^ 

~ 2^^ ' ^ ^ 
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9± = 9±± TT : l=rr, 2 M2 I /I ^/Q213/2 ^ ~ > [^^) 



1 

7± = 



2 



A, + A2 ± v/^(l + a^/.)-(g + .) _ 



(21) 



Note that these equations can be obtained from Eqs. (O U\ [El [El and [TT]) with the substi- 
tution (g ± i/a) — [g + k ± ^/a{l + a^n)]. Additional sources of precession, such as stellar 
and planetary equatorial bulges, will add additional terms to the diagonal of the matrix A 
and can be treated in a similar manner. 

The relativistic effect increases both secular rates (Eq. [T8|) since it causes the orbits 
to precess in the same direction as the secular interaction does. The mode eccentricity ratio 
|?7^| (Eq. [T9|) increases significantly with GR, or with any source of precession that favors the 
inner orbit, indicating that it is more difficult to force the eccentricity of a rapidly precessing 
inner orbit in the aligned mode. The ratio |?7f_|, however, decreases slightly; increasing the 
precession of the inner orbit in the anti-aligned mode actually strengthens secular coupling. 
As for the mode damping rates (Eq. fIT\i . relativistic precession decreases the aligned mode 
damping rate, but increases that of the anti-aligned mode. It also decreases e, the deviation 
angle of the mode apsidal lines from perfect alignment (Eq. 



3. Applications to the Observed Systems 

We now apply the theory developed above to close-in exoplanets, beginning with known 
two-planet systems. We then proceed to systems in which there is linear trend in the star's 
radial velocity that may signal the presence of an outer companion and finally we consider 
systems with no hint of a companion. 



Table 1: Properties of planets and stars discussed in this paper 



Planet 


nip {mj) 


Rp (Rj) 


a (AU) 


e 


con 


m* (m0) 


Age (Gyrs) 


Te (Gyrs) 


HAT-P-13b 


0.85 ±0.0354 


1.28 ±0.079 


0.0427 ±0.000875 


0.013 ±0.0041 




-1 9+0.05 
-■^•^-0.1 


c:+2.5 
^-0.8 


- 0.125 


HAT-F-lSc 


14.3 ±0.691 




1.23 ±0.0251 


0.662 ± 0.0054 


175.3 ±0.35 


1 9+0.05 
^■^-0.1 


c;+2.5 
"J-0.8 




HD 187123b 


0.51 ±0.0173 


1* 


0.0421 ±0.000702 


0.010 ±0.00593 


24.5 


^•^^-0.024 


5.33 


- 0.3 


HD 187123 c 


1.9 ±0.152 




4.8 ±0.367 


0.25 ±0.0334 


240 ± 18.6 


1 04+0 026 
^•"^-0.024 


5.33 




GJ436b 


0.073 ±0.00318 


n 07^7+0.0082 

U.O(U(_g 0092 


0.0287 ±0.000479 


0.16 ±0.019 


351 ± 1.2 


.r+0.014 

^'•'±0-0.012 




- 2 


BD-103166b 


0.43 ±0.0174 


1* 


0.0438 ±0.000730 


o.o21e;°42 


334 


n Q9+0.046 
'-'•-'^-0.024 


4.18 


- 0.42 


HAT-P-26 b 


0.059 ±0.00718 


0.565 ±0.052 


0.0479 ±0.000798 


0.12 ±0.06 


100 ± 165 


0.82 ±0.033 


q+3 


- 2.6 


WASP-34b 


0.58 ±0.0285 


1 99+0.11 
^■^^-0.08 


0.052 ±0.00120 


0.04 ±0.0012 


320 ±20.9 


1.01 ±0.07 


6 7+6-^ 


- 0.7 


HD 1491 43 b 


1 33 ± 0784 


1* 


053 ± 00147 


016 ± 01 





~ 1.2 ± 0.1 


7.6 ± 1.2 


~ 2.75 


HAT-P-21b 


4.1 ±0.173 


1.024 ±0.092 


0.0495 ±0.000825 


0.23 ±0.016 


309 ±3 


0.95 ±0.042 


10.2 ±2.5 


- 6 


HAT-P-23 b 


2.1 ±0.122 


1.368 ±0.09 


0.0232 ±0.000387 


0.11 ±0.044 


120 ±25 


1.13 ±0.035 


4± 1 


- 0.005 


HAT-P-32 b 


1.0 ±0.169 


2.037 ±0.099 


0.0344 ±0.000574 


0.16 ±0.061 


50 ±29 


-1 -1 O+0.043 


q + 1.5 


- 0.005 


HAT-P-33 b 


0.8± 0.117 


1.83 ±0.29 


0.050 ±0.00115 


0.15 ±0.081 


100 ± 119 


1.40 ± 0.096 


2.4 ±0.4 


- 0.06 


HD 88133 b 


0.30 ± 0.0270 


1* 


0.0472 ±0.000786 


0.13 ±0.072 


349 


1.2 ±0 


9.56 


- 0.38 
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3.1. Two-Planet Systems 



There are more than 30 multi-planet systems that host one or more close-in planets 
with the orbital period Porb ^ 20 days. Out of these, many systems including Gliese 876, 55 
Cnc, and v And have three or mo re planets, which makes apsidal analysis more complicated 
(e.g., iBarnes fc Greenberg 



2006al ). Moreover, most Kepler-detected planets have unknown 
eccentricities, and are not the best candidates for our analysis. After removing multiple 
planet systems and those with poorly determined eccentricities, we are left with eight two- 
planet systems, with which we can test our theory. 

First, we compute eccentricity damping timescales for the inner planets (rg) by mak- 
ing assumptions about the strength of the tidal interactions. We then compare ije) to the 
stellar age (rAgc) to determine whether an orbit has had sufficient time to circularize. For 



all systems, we adopt conventio nal planetary and stellar tidal quality 



= 10^, respectively (see e.g.. 



Jackson et al. 



20081: 



Matsumura et al. 



'actor s = 10^ and 



2010 



and references 



therein). We convert these to modified tidal quality factors with Q' = 1.5Q/k2. For giant 
planets in our Solar System (Jupiter, Saturn, Uranus, and Neptune), the measured gravi- 



tational moments agree well with an n ~ 1 polvtrop e (IHubbard 



which corresponds to the Love number k2 = 0.52 ( iMotz 



other hand, then =_3_polytrope is a good approximation (e.g.. 



k2 = 0.028 flMotz 



1974 



Bobrov et al. 



19781 ) 



1952). For most stars, on the 



Horedt 



2004 - which yields 



1952h . We finally obtain Q; = 2.88 x 10^ and Q', = 5.36 x 10^. We 
estimate the eccentricity damping timescale Tg by integrating a set of tidal equations based 
on the equilibrium tide model from the measured planetary eccentricity down to e = 10~^, 
and further assume that the tid al quality factors evolve proportional to the inverse of the 



mean motion Q oc l/rii (see e.g.. 



Matsumura et al. 



20101). Some of our results appear in the 



final column of Table [T] Note that these simulation results compare favorably to the simpler 
approximation of Eq. [T^when one properly accounts for the multiple e- folding times needed 
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to damp eccentricities to e = 10~^. Only 2 of the 8 two-planet systems (HD 187123 and 
HAT-P-13) have a short tidal circularization times Tg < rAgc; physical and orbital parame- 
ters of these two systems can be found at the top of Table [TJ Since « TAge, the systems 
are likely to have been significantly modified by tides. 

Next, we check the strength of secular interactions for the two systems. In the HD 187123 
system, the inner planet is at ~ 0.04 AU while the outer one is at ~ 5AU. With this 
configuration, the inner planet is strongly bound to the star and its interaction with the 
outer planet is weak. The inner and outer planets of HAT-P-13, however, are more closely 
spaced at ~ 0.04 and ~ 1.2 AU, respectively. Furthermore, since m2 >> mi and 62 = 0.66 
is large (Table [1]), the secular forcing of the inner planet is substantial. Thus, out of the 
eight systems, only HAT-P-13 has both strong secular interactions and a short tidal damping 
time; it is, accordingly, the best test case for our theory. 

We now utilize the linear secular theory without tides developed in Section [2?n and show 
that the model predicts the current apsidal state of HAT-P-13 reasonably well. In Fig. [TOl we 
compare the apsidal state of HAT-P-13 estimated by sec ular theory with that obta ined from 



2002). Despite 



a direct N-body simulation done with the HNBody code (iRauch fc Hamilton 
that HAT-P-13c's large eccentricity of ~ 0.66 violates the assumption of our basic secular 
theory, discrepancy between the two integrations is merely ~ 25% (Fig. [TOj) . It appears that 
the system librates with large amplitude about Aw = 0° so that the system is not far from 
the aligned separatrix (Fig. 12.31) . This is yet another example of a ge neral observation that 



extrasolar multi-planet systems tend to be near a secular separatrix (IBarnes fc Greenberg 



2006bt). The near-separatrix state of this system is somewhat surprising given that the 
tidal decay time is a tiny fraction of the stellar age (Table [1]) and so the system should 
have long ago damped to the apsidally-locked state. Perhaps, given the uncertainty in Q', 
our estimation of Te is off by a large factor. Alternatively, perhaps there is a third as-yet 



- 30 - 




Fig. 10. — The current apsidal state of the HAT-P-13 system. We show the orbits of both 
planets on the e exp{iAw) plane, similar to the lower plots in Fig. [7] and [91 The orange curves 
represent the solution of secular equations, while blue curves are obtained from an N-body 
simulation. The orbits librate about Aw = 0°, with a libration amplitude ~ 67° predicted 
by secular theory and ~ 50° measured from the N-body simulation. The narrowness of the 
outer planet's arc is both cases is due to the fact that m2 >> mi. 
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undiscovered planet affecting this secular system. 

We next compare the current state of the system with the expectation from tidal- 
damping theory as developed in Section 12.21 We assume that the parameters of the inner 
planet's orbit are known, assume further that one mode has fully damped away (despite 
the contrary evidence of Fig. [TT]) . and proceed to predict parameters of the outer planet. 
We begin by asking which mode is favored which depends on the damping rates given by 
Eq. (fT6|) . or Eq. (12T]) when GR effects are important. 

When a system is locked into one of the secular eigen-modes, the outer and inner orbits 
have a predictable eccentricity ratio |?7-|-| = {e2/ei)± (see Eqs. [T71 and [221) . In Fig. [TTl we 
show contour plots of the predicted eccentricity of the outer planet 62 for the more slowly 
damped eigen-mode in the parameter space (g, a). The dashed curve divides the space into a 
region in which the slow anti-aligned mode persists (top left, red area) and a region in which 
the slow aligned mode survives (bottom right, blue area). In the shaded areas, the lifetime 
of the slower mode is longer than the age of the system. Conversely, in the white central 
area, both modes should have already damped away and both orbits would be circular by 
now. Given the inner planet's non-zero ei = 0.013 (Tabled]) and the fact that Tg << TAge, 
we expect the outer planet to be in one of the shaded regions. Furthermore, since e > 1 
correspond to unbound orbits, these parts of the shaded areas in Fig. [11] are also off limits. 

The boundaries of these colored areas are determined by equating the age of the system 
(5 Gyr for IIAT-P-13) to the circularization time Tg in Table [T] An older system age TAge 
and/or faster damping timescale Te would expand the white area outward away from the 
dashed line. Conversely, if tides have not been active over the full stellar age, as is possible 
for a recent resonance crossing or a planet-planet scattering event, the white area would 
shrink inward toward the dashed line. 
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q = tTij/m-i 



100.00 



0.01 0.10 1.00 10.00 100.00 

q = tTij/m, 



a) without general relativity b) with general relativity 

Fig. 11. — Orbital states of possible companions for HAT-P-13b without (a) and with (b) 
general relativistic precession. Each point in these q - a plots represents an outer companion 
with corresponding mass and semi-major axis. The dashed curve divides the plane into 
regions in which the aligned mode damps faster (top-left) and the anti-aligned mode damps 
faster (bottom-right). Furthermore, in the shaded regions, either the aligned mode (blue 
shading) or the anti-aligned mode (red shading) can survive tidal dissipation and last longer 
than the age of the system (5 Gyrs). Finally, the solid contour lines represent the eccentricity 
of the outer planet assuming that it is in the long-lived mode. The location of HAT-P-13c 
is marked by a white symbol with error bars. 
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The effect of the relativistic pericenter precession is substantial and can be quantified 
by comparing Fig. [TTk and b. Relativistic perturbations are large for close-in planets and 
can dominate planetary perturbations when the outer companion is less massive (small q) 
and/or distant (small a). As discussed in section [531 adding GR diminishes the anti-aligned 
area (red shading) due to faster damping rates but significantly expands the aligned area 
(blue shading) in accordance with Eq. fl^T]) . Notice that, as expected, the lower left quadrant 
of the plot, small q and small a, experiences the greatest changes. Conversely, the changes 
to the shading of the other three quadrants is relatively minor. The dashed lines in Fig. [TTb 
and b are given by g = a/o and q + k= y/oi{l + o^k), respectively; these expressions simply 
compare the real diagonal elements of the respective A matricies. The anti-aligned mode 
damps most quickly if the inner planet precesses faster while the aligned mode damps first 
if the outer planet precesses faster. Along this dashed line, the difference between the mode 
precession rates ((^^ — gt) is minimized and the mode damping rates 7± = (Ai + A2)/2 are 
identical. 

Finally, note that adding GR significantly changes the contours for the outer planet's 
eccentricity. In the region of aligned libration, the eccentricity contours are moved upward 
by the inclusion of GR, indicating that a more eccentric outer planet is required to maintain 
the mode, for a given a and q. Conversely, in the anti-aligned region, the contours move 
to the left indicating that a lower eccentricity on the outer planet is needed to preserve the 
mode. The reasons for these changes were discussed in section 12.41 

The actual location of the outer planet HAT-P-13c is marked in both panels with white 
error bars; notice first that the different predictions for 62 indicate that GR cannot be 
neglected for this system. The planet resides well within the more slowly damped aligned 
mode region, as expected from linear secular theory (see also Fig. [10]). Figure [TTb suggests 
that the eccentricity of the outer planet is close to 1, while the observed value is actually 
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0.662 ±0.054 (Tabled]). Although the agreement is not perfect, our method does successfully 
predict apsidal alignment and a large eccentricity for the outer planet. In fact, our theory, 
derived for small eccentricities only, does surprisingly well outside of its realm of strict 
validity. We conclude that the secular perturbations between the two planets in the HAT- 
P-13 system are most likely responsible for the non-zero eccentricity of the inner planet. 

In the following two sections, we study single-planet systems with and without linear 
trends in the stellar radical velocity, and investigate whether their non-zero eccentricities 
could be explained by the potential companions. 

3.2. Single-Planet Systems with a Hint of a Companion 

As shown in Section [272] and [273| a two-planet system should have evolved into either an 
aligned or an anti-aligned apsidally locked state when the tidal dissipation is strong enough. 
Equations flTI?]) and (^T^ thus provide a single constraint on the three parameters of the 
unknown outer companion: the mass ratio q, the semi-major axis ratio a, and the eccentricity 
ratio of the two planets. Therefore, we can predict a range of possible companions that may 
force a non-zero eccentricity on an observed close-in planet. We illustrate our method with 
several examples here and in the next section. 

There are 16 single, close-in planet systems with Porb < 20 days and non-zero eccen- 
tricities that have an observed linear trend in the stellar radial velocity, which indicates 
the possible existence of a companion on a more distant orbit. For these systems, we can 
place a unique constraint on the potential companion. For simplicity, we exclude the two 
systems that have a large projected stellar obliquity (HAT-P-11 and WASP-8), and compare 
the estimated Tg with the stellar age TAgc for each of the remaining systems, as described 
above. We find that five of fourteen systems (GJ436, BD -10 3166, HAT-P-26, WASP-34, 



- 35 - 



and HD 149143) have Tg < TAgc- We list parameters for these systems in the middle section 
of Table [1] As the properties of the putative companions are unknown, we cannot test for a 
strong secular interaction as in the previous section and so we investigate all five systems. 

Figure [T2k is similar to Figure [TTb. but for the planetary system around GJ436. The 
thick dashed curve is an upper limit to the outer planet's mass estimated from the observed 
linear trend. Here, we simply assume that the minimum mass of a potential outer planet 
i s expressed as m2 = a\ait/G, where an = 1.36 ± 0.4ms~^yr~^ is an observed linear trend 



( iManess et al. 



2007l )lH The dotted curve indicates an observation limit for the radial velocity 
(RV) method. The solid portion of this curve is plotted as a reference, and shows the limit 
estimated for the 1-year of the observation period. To the right of this solid curve, a full 
orbit of a hypothetical outer planet is observable within a year. To plot this, we express the 
mass of a potential outer planet as m2sin2 = m ^,v^, sin i / ^yGm^Ja^, where i is the viewing 
angle, and assume the RV limit of v^, sini = lms~^. 

The observed planet GJ436b is about a Neptune-mass object (msinz = 0.073 mj ~ 
1.35 mjv) which is 0.0287 AU away from the central star (orbital period ~ 2.6 days), and has 
an orbital eccentricity of 0.16 ± 0.019 (Tabled]). If this eccentricity is due to another planet 
and the system has damped to an eigen-mode, then eccentricity contours in Fig. [T2k show 
that broad aligned and anti-aligned regions are allowed for a potential companion, except for 
small (m2 < 0.2 m^v) and/or distant (a2 > 0.24 AU) planets. Furthermore, the plotted RV 
limit indicates that nearly all hypothetical outer planets which could be responsible for the 
high eccentricity of G J 436 b should be observable within a year. The curve representing the 
maximum linear trend, however is far below the e = 1 contour, implying that this potential 
planet cannot be responsible for the current eccentricity of the observed planet; given its 



^This long-term trend has not been confirmed by HARPS (jBonfils et al 



20131) 
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Fig. 12. — Eccentricities of possible companions for Hot Jupiters a) GJ436, and b) WASP-34 
including GR effects as in Fig. ITTb . In the white area, both modes should damp away within 
the age of the system of 6 and 6.7 Gyr, respectively. The central stars of these systems each 
have an observed linear trend in their radial velocities, 1.36 and 55 ms~^yr~^ respectively, 
that may be indicative of second planets. These constraints are given by the thick dashed or 
solid-dashed curves. For WASP-34, the region to th e right of the thick so lid curve represents 



201l|). The solid-dotted 



where a potential candidate is expected to exist (iSmalley et al. 
curves represent the RV observation limit of 1 ms~^; only planets to the right of this curve 
are detectable with current technology in a 1-year observational period. 
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great distant, its secular interactions are simply too weak. This re sult is consistent with th e 
comparison of the secular and tidal circularization timescales by 



Matsumura et al. 



(120081). 



Since the system has a tidal dissipation timescale {t^ ~ 2 Gyr) comparable to the stellar age 
(^Agc ~ SlgGyr), a non-zero eccentricity of this planet might also be explained within the 
uncertainties of the stellar age. 

Another example is shown in Figure [T2b for WASP-34, which has an observed planet 
of msini = 0.58mj, a = 0.052 AU, and e = 0.038 ± 0.012 (see Tabled. Again, most of 
the parameter space of the q — a plane is available for a possible secular companion, except 
for small (m2 < 9.2 m^;) and/or distant (02 ^ 1.3 AU) planets. The solid portion of the RV 
limit again indicates that such a companion should be observ able within a year in most cases. 



20111 ). Since the 



The system has an observed linear trend of 55 ± 4 m s^^yr"^ (iSmalley et al.l 
long-term trend has not reached its maxima or minima, the orbital period of the outer body 
has to be grea ter than twice the R V data baseline. This corresponds to 02 ^ 1.2 AU and 



m2 > 0.45 mj (jSmalley et al. 



2OIII ). and is indicated as a thick solid curve in the figure. The 



linear trend lies much closer to the critical e = 1 contour than GJ436. The uncertainties 
in both observations and the secular model could bring these curves closer. A companion 
planet responsible for both this trend and the eccentricity of the WASP-34 b will likely have 
a very high eccentricity. 



If no companion is present, how do we explain WASP-34? For this system, t 



le eccentric 



20111 ) and. 



orbit model gives only a slightly better fit than the circular one (ISmalley et al. 
since the estimated eccentricity damping time (rg ~ 700 Myr) is short compared to the 
stellar age (rAgc ~ ^-^-t'lGyi), the inner orbit might well be circular. Alternatively, if the 
orbit is truly eccentric, we would need to assume about an order of magnitude less efficient 
tidal dissipation to explain this system. Finally, the system could also have undergone some 
dynamical event lately which changed the original eccentricities. 
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The other systems with a hnear trend (BD -10 3166, HD 149143, and HAT-P-26) show 
a similar result to GJ436 (Fig. [T2b). and thus a potential planet is too far to force the 
eccentricity of the inner planet to its current value. Interestingly, the observed eccentric- 
ities are low and consistent with zero for BD -10 3166 and HD 149143 ( 



Fischer et al 



20061 ). and poorly constrained for HAT-P-26 flHartman et al 



Butler et al 



2000 



2011a). BD- 



10 3166 has a short tidal dissipation timescale (rg ~ 420 Myr) compared to the stellar age 
(TAgc ~ 4.18 Gyr), so the circular orbit assumption makes sense. On the other hand, both 
HD 149143 and HAT-P-26 have relatively long dissipation timescales (rg ~ 2.75 Gyr and 
~ 2.6 Gyr, respectively) compared to the stellar ages (rAgc ~ 7.6 ± 1.2 Gyr and ~ 9;'^4 g Gyr, 
respectively). Thus uncertainties in the age estimates and/or in the tidal dissipation rates 
could allow close-in planets to maintain their eccentricities without assistance. 

In summary, we have found no single-planet systems, where the none-zero eccentricities 
could be explained by perturbations from hypothetical planets corresponding to the observed 
linear trends. 



3.3. Single-Planet Systems with no Hint of a Companion 

Since there are many single, close-in planet systems without linear trends, we focus 
on planets whose non-zero eccentricities are hardest to explain — the closest-in exoplanets. 
There are 8 single planet systems with an orbital period Porb < 5 days, eccentricity > 0.1, 
and a small or unknown stellar obliquity. Comparing the tidal timescale to the stellar age, we 
find that five of eight such systems have < ta^c (HAT-P-21, HAT-P-23, HAT-P-32, HAT- 
P-33, and HD 88133, see Tabled]). All of these are hot Jupiter systems. For the remaining 
three systems, the stellar ages of KOI-254 and Kepler-15 are unknown, while GJ674 has a 
very long Tg > 10 Gyr compared to the stellar age (~ 0.55 Gyr). 
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Fig. 13. — Eccentricities of possible companions for Hot Jupiters a) HD 88133, and b) HAT- 
P-21 including GR effects as in Fig. [TTb . In the white area, both modes should damp away 
within the age of the system of 9.56 and 10.2 ± 2.5 Gyr, respectively. As the inner HAT-P- 
21 planet has a long damping time of about 6 Gyr, all eccentricities can be sustained and 
there is almost no white area visible. The thick dashed curves represent the RV observation 
limit of lms~^, with the solid portion indicating a 1-year observing period. All potential 
HD 88133 outer planets reside in the aligned zone are detectable. Both aligned and anti- 
aligned solutions exist for HAT-P-21 and nearly all outer planets are detectable. 
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Figure [T3k is similar to Figure [TTb . but for the planetary system around HD 88133. 
The planet is 0.0472 AU away from the central star with an orbital period of 3.4 days and 
an eccentricity of 0.13 ± 0.072 (see Tabled]). Since the circularization time is estimated to 
be very short (rg ~ 380 Myr) compared to the stellar age (TAgc ~ 9.56 Gyr), a moderately 
high eccentricity of this planet is surprising. Figure [T^ excludes the entire anti-aligned 
region for a potential companion. The figure also excludes most companions with small 
mass and/or long orbital period — more precisely, those with m2 < 0.2 mi ~ 0.06 mj 
and/or 02 > lOai ~ 0.472 AU. This is understandable since these planets would have weak 
secular interactions with the inner planet. Thus, a potential companion is expected to be in 
the aligned region, massive, and close to the star. As we can see from the RV observation 
limit, such massive planets in the aligned libration region would be observable within a year, 
although no such planet has been found. HD 88133 has a low stellar jitter (~ 3 .2 ms"^), and 



the e ccentric orbit assumption works only slightly better than the circular one (IFischer et al. 



20051 ). Thus, unless the tidal quality factor for this system is very different from what we 
have assumed here, we argue that the true eccentricity of HD 88133 b is actually near zero. 
Follow-up observations would better constrain the eccentricity and the existence or absence 
of a potential companion for this system. 

The case for HAT-P-21 is shown in Figure [T3b . The planet is 0.0495 AU away from 
the central star (orbital period ~ 4.1 days), and has an orbital eccentricity of 0.23 ± 0.016 
(Tabled]). Figure dSb allows a very broad parameter space for a possible secular companion; 
in fact, the white zone indicating efficient eccentricity damping is nearly absent. But this 
undercuts our assumption that the system has had time to damp into a pure eigenmode 
and, accordingly, the eccentricity contours are not reliable. If we proceed with the dubious 
assumption of a single mode, the figure does not allow most companions with m2 < 0.82 mj 
and/or 02 ^ 0.495 AU. Furthermore, the Im/s radial velocity limit observationally precludes 
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almost any planet that can be significantly coupled to HAT-P-21 b. Accordingly, we seek 
another explanation for the eccentricity of this system; it can be naturally explained if tidal 
dissipation were just slightly less efficient than we have assumed, since the estimated tidal 
circularization time is relatively long [r^ ~ 2.6 Gyr) compared to its stellar age (TAgc ~ 
QtigGyr). 

The other systems (HAT-P-23, HAT-P-32, and HAT-P-33) have a similar trend to 
HD 88133 (Figure [T3h). with the entire anti-aligned region being excluded for a potential 
companion. Since all of their circularization times are more than 1-2 orders of magnitude 
shorter than the estimated stellar ages, the moderately- high eccentricities (e > 0.1) of these 
planets need to be explained. From figures similar to Fig. [131 "we find that potential com- 
panions for these systems tend to be more massive than the observed planets, and thus are 
likely to be observable by the radial velocity method. However, no companio ns have been 



found. It is interesting to no te that all of these systems have high stellar jitters (IBakos et al. 



2011 



Hartman et al. 



2011bl ). Although this could mean that potential companion planets 
are difficult to observe, high jitters also lead to poorly-constrained orbital eccentricities. In 
fact, all of these planets can also be fit well with the circular orbit model. Our model suggests 
that the circular orbits are probably the most likely solution. Future observations that yield 
a more accurate solution for the eccentricity are sorely needed. 

For our analysis of exoplanetary systems in this section, we did not explicitly take 
account of the effects of uncertainties in orbital or stellar parameters. In particular, errors 
in stellar ages and eccentricities are often large and may change our results significantly. 
We have tested such effects for all the systems we discussed in this section, and found that 
our conclusions will not change within the currently estimated uncertainties in parameters. 
Also, the assumption of an apsidal lock that we made in Section 13.21 and 13.31 may be too 
strong; instead it is possible that not enough time for complete damping of one secular 
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mode has elapsed. In this case, for a close-in planet with known mass, semi-major axis, and 
eccentricity, the constraint on 1712, 02, and 62 is rather approximate than exact. 

4. Discussions and Conclusions 

The eccentric orbits of single, close-in planets are generally circularized on timescales 
shorter than the stellar ages. On the other hand, close-in planets in multiple-planet systems 
have much longer tidal circularization timescales and thus are able to maintain eccentric 
orbits for the stellar ages or longer. Given this difference in tidal circularization times, we 
might expect a difference in the eccentricity distributions of close-in planets with and without 
known companions. The eccentricities forced by secular perturbations from an outer planet 
are typically quite small, however, and so it is perhaps not surprising that no such difference 
has yet been observed. 

In this paper, we have explored the possibility that the non-zero eccentricities of close-in 
planets are due to observed or hypothetical planetary companions. We have provided an 
intuitive interpretation of a simple secular evolution model of a coplanar two-planet system 
that includes both the effect of the orbital circularization (Section 12. 2p and GR corrections 
(Section 12.41) . We have tested our model by comparing the evolution of apsidal states and 
orbital eccentricities with N-body simulations, and found that the agreement between the 
model and the simulations is very good. We have also applied our model to all of the relevant 
two-planet systems (Section 13.11) . as well as single-planet systems with and without a long- 
term trend in the radial velocity to indicate a possible second planet (Sections 13.21 and 13. 3[ 
respectively). Our main results are summarized below. 

1. In the lowest-order secular theory, the evolution of non-dissipative two-planet systems 
is described by a linear combination of two modes characterized by pericenter alignment 
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and ant i- alignment. Eccentricity damping slightly shifts the two normal modes from 
perfect symmetry, which speeds up the precession rate of the aligned mode and slows 
that of the anti-aligned mode (see Section [272]) . 

2. Eccentricity damping affects the two modes at different rates. Accordingly, the ap- 
sidal state of a two-planet system transitions between libration and circulation, and 
eventually is locked to either an aligned or anti-aligned state (see Section [2.3p . The 
eccentricity of both planets subsequently decays at a very slow rate. 

3. General relativistic effects increase the precession rates of both aligned and anti-aligned 
modes. As a result, they decrease the aligned-mode damping rate, and increase the 
anti-aligned mode damping rate (see Section [22 



We have confirmed results o 



Greenberg fc Van Laerhoven 



previous studies (e.g. 



Wu fc Goldreich 



2002 



Mardling 



2007 



201 ll ). and shown that close-in planets in multiple-planet sys- 
tems can maintain non-zero orbital eccentricities substantially longer than the single ones. 
We found, however, that there are currently few two- or one-planet systems which show signs 
of secular interactions that are strong enough to significantly slow tidal circularization. In 
Section [STT] we found that one out of eight systems (HAT-P-13) shows evidence that secular 
interactions are slowing orbital circularization. In Section 13. 2^ we apply our model to 14 
single-planet systems with a linear trend, and found that five of fourteen have Tg < TAge- 
Our secular model predicts that none of their eccentricities are likely to be affected by hy- 
pothetical planets which could cause the long-term linear trends. We have further studied 
eight very close-in (Porb < 5 days), significantly eccentric (e > 0.1) single-planet systems in 
Section 13.31 We found five of eight systems which could not be explained by a single-planet 
orbital circularization with the conventional tidal quality factors. Potential companions for 
all of these systems are massive planets in the apsidally-aligned region, and should be observ- 
able with current technology. Since all of the host stars have high stellar jitters, it is possible 
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that the planetary eccentricities are systematically overestimated or that outer planets are 
more difficult to observe than we have assumed here. 

There are some limitations and caveats in our model. We have adopted the leading- 
order secular theory for two-planet systems, and in principle our model cannot be applied 
to high eccentricity or high inclination systems. However, the model predicts the general 
trend of apsidal states fairly well even for a highly eccentric case (see Section ISTTj) . More- 
over, recent observa. t ions indicate that mu ltiple-planet systems tend to be well-aligned (e.g.. 



Figueira et al. 



2012 



Fabrycky et al. 



2OI2I ). Nevertheless, it would be useful to extend this 
kind of a study to higher eccentricities and inclinations, and also to systems with more than 
2 planets. Also, we have ignored the slow decay in semimajor axis due to tides, which should 
be weaker than eccentricity damping by a factor of e^; this is consistent with the low e 
assumption made by linear secular theory. Nevertheless, previous studies have shown that 



inward migration (e.g.. 



Wu & Goldreich 


2002; 


Greenberg & Van Laerhoven 


2011) 



consequence of different damping rates is that the eccentricity ratio does not remain constant 
as the apsidally-locked state evolves. Our expression f or the mode damping rates (Eq. |2T|) is 



consistent with that of 



Greenberg fc Van Laerhoven! (120111 ) (see their Equation 19), in the 



limit of no migration. The tidal and rotational deformations of planets and stars also change 
orbital precession rates, and these effects may become more important than the GR effect 
for very close-in planets or rapidly spinning st ars. These effects c an be easily added to our 
model using the techniques of section [531 (e.g.. 



Laskar et al 



20121 ) 



Overall, our study indicates that secular interactions slow down the tidal circularization 
of the inner planet while speeding up that of the outer planet. Our survey of likely systems 
available to us in 2012 indicate that secular interactions may not be a dominant cause for 
the currently observed hot, eccentric planets. The lack of close-in planets with strong secular 
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interactions may be par tially explained by inward orbital decay that is accelerated by non- 



zero eccentricities (e.g. 



Adams fc Laughlin 



2006bl ). It remains for future observations to 



determine whether the scarcity of compact secular systems will persist. With improved 
statistics and precision measurements of close-in exoplanet eccentricities, it may become 
possible to find diagnostic differences in the eccentricity distributions for single- and multiple- 
planet systems. 



REFERENCES 

Adams, F. C, & Laughlin, G. 2006a, ApJ, 649, 992, 992 
— . 2006b, ApJ, 649, 1004, 1004 

Bakos, G. A., Hartman, J., Torres, G., et al. 2011, ApJ, 742, 116, 116 
Barnes, R., & Greenberg, R. 2006a, ApJ, 652, L53, L53 
— . 2006b, ApJ, 638, 478, 478 

Batygin, K., Laughlin, G., Meschiari, S., et al. 2009, ApJ, 699, 23, 23 
Batygin, K., & Morbidelli, A. 2013, A J, 145, 1, 1 

Beauge, C., Ferraz-Mello, S., & Michtchenko, T. A. 2003, ApJ, 593, 1124, 1124 

Bobrov, A. M., Vasil'Ev, R R, Zharkov, V. N., & Trubitsyn, V. R 1978, Soviet Ast., 22, 
489, 489 

Bonfils, X., Delfosse, X., Udry, S., et al. 2013, A&A, 549, A109, A109 
Brouwer, D., van Woerkom, A., & Jasper, J. 1950 

Burns, J. A. 1977, in Planetary Satellites, ed. J. A. Burns (Tuscon, AZ, USA: Univ. of 
Arizona Press), 113-156 

Butler, R. P., Marcy, G. W., Fischer, D. A., et al. 1999, ApJ, 526, 916, 916 

Butler, R. P., Vogt, S. S., Marcy, G. W., et al. 2000, ApJ, 545, 504, 504 

Butler, R. P., Wright, J. T., Marcy, G. W., et al. 2006, ApJ, 646, 505, 505 



- 46 - 



Chiang, E. 1. 2003, ApJ, 584, 465, 465 
Danby, J. M. A. 1988 

Fabrycky, D. C, Lissauer, J. J., Ragozzine, D., et al. 2012, ArXiv e-prints, arXiv:1202.6328 
Figueira, P., Marmier, M., Boue, G., et al. 2012, A&A, 541, A139, A139 
Fischer, D. A., Laughhn, G., Butler, P., et al. 2005, ApJ, 620, 481, 481 
Fischer, D. A., Laughlin, G., Marcy, G. W., et al. 2006, ApJ, 637, 1094, 1094 
Goldreich, R. 1963, MNRAS, 126, 257, 257 

Greenberg, R. 1977, in Planetary Satellites, ed. J. A. Burns (Tuscon, AZ, USA: Univ. of 
Arizona Press), 157-168 

Greenberg, R., & Van Laerhoven, C. 2011, ApJ, 733, 8, 8 

Hamilton, D. P. 1994, Icarus, 109, 221, 221 

Hartman, J. D., Bakos, G. A., Kipping, D. M., et al. 2011a, ApJ, 728, 138, 138 

Hartman, J. D., Bakos, G. A., Torres, G., et al. 2011b, ApJ, 742, 59, 59 

Horedt, G. P., ed. 2004, Astrophysics and Space Science Library, Vol. 306, Polytropes - 
Applications in Astrophysics and Related Fields 

Hubbard, W. B. 1974, Icarus, 23, 42, 42 

Jackson, B., Greenberg, R., & Barnes, R. 2008, ApJ, 678, 1396, 1396 
Ketchum, J. A., Adams, F. C., & Bloch, A. M. 2013, ApJ, 762, 71, 71 
Laskar, J., Boue, G., & Correia, A. C. M. 2012, A&A, 538, A105, A105 
Lee, M. H. 2004, ApJ, 611, 517, 517 

Maness, H. L., Marcy, G. W., Ford, E. B., et al. 2007, PASP, 119, 90, 90 
Mardling, R. A. 2007, MNRAS, 382, 1768, 1768 

Matsumura, S., Peale, S. J., & Rasio, F. A. 2010, ApJ, 725, 1995, 1995 
Matsumura, S., Takeda, G., & Rasio, F. A. 2008, ApJ, 686, L29, L29 
Motz, L. 1952, ApJ, 115, 562, 562 



Murray, C. D., & Dermott, S. F. 1999 

Namouni, F. 2007, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 729, Lecture 
Notes in Physics, Berlin Springer Verlag, ed. D. Benest, C. Froeschle, & E. Lega, 233 

Peale, S. J. 1986, in Satellites, ed. J. A. Burns & M. S. Matthews (Tuscon, AZ, USA: Univ. 
of Arizona Press), 159-223 

Pont, F., Husnoo, N., Mazeh, T., & Fabrycky, D. 2011, MNRAS, 414, 1278, 1278 

Rasio, F. A., Tout, C. A., Lubow, S. H., & Livio, M. 1996, ApJ, 470, 1187, 1187 

Ranch, K. P., & Hamilton, D. P. 2002, in Bulletin of the American Astronomical Society, 
Vol. 34, Bull. Am. Astron. Soc, 938 

Shen, Y., & Turner, E. L. 2008, ApJ, 685, 553, 553 

Smalley, B., Anderson, D. R., Collier Cameron, A., et al. 2011, A&A, 526, A130, A130 

Winn, J. N., Johnson, J. A., Howard, A. W., et al. 2010, ApJ, 718, 575, 575 

Wright, J. T., Upadhyay, S., Marcy, G. W., et al. 2009, ApJ, 693, 1084, 1084 

Wu, Y., & Goldreich, P. 2002, ApJ, 564, 1024, 1024 

Zakamska, N. L., Pan, M., & Ford, E. B. 2011, MNRAS, 410, 1895, 1895 

Zhang, K., & Hamilton, D. P. 2007, Icarus, 188, 386, 386 

— . 2008, Icarus, 193, 267, 267 



This preprint was prepared with the AAS IATJt;X macros v5.2. 



